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We describe a practical procedure to calculate the Coulomb matrix elements of 2D spatially 
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important condensed matter finite systems. We derive an analytical expression, for arbitrary sepa- 
rations, in terms of a single infinite series and apply a u-type Levin transform in order to accelerate 
the resulting infinite series. This procedure has proven to be efficient and accurate. Direct conse- 
quences concerning the functional dependence of the matrix elements on the separation distance, 
transition amplitudes and the diagonalization of a single electron-hole pair in vertically stacked 
parabolic quantum dots are presented. 
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The interesting physical properties displayed by charge 
carriers in spatially confining settings promise to have a 
diverse range of applications relevant to technological de- 
velopments in the near future. Systems such as semicon- 
ductor quantum wells, quantum wires and quantum dots, 
will have great impact in the progress of quantum compu- 
tation and quantum information [1] and in producing new 
technological devices in electronics, spintronicsQ and op- 
toelectronics [3]. Furthermore, on the experimental and 
theoretical side, systems which have acquired increasing 
relevance over the past decade, such as indirect excitons 
B, graphene bilayers p] and indirect magnetobiexcitons 
p , show great promise for present-day and future under- 
standing of condensed matter systems. 

In order to grasp the full potential of these advances, it 
is necessary to have an increasingly detailed theoretical 
understanding of the physics of confined charge carriers. 
These studies have included, for example, learning about 
their optical properties, describing their response to ex- 
ternally applied electric fields, and understanding their 
behaviour when they are subject to external magnetic 
fields 0,0. 

Such theoretical studies have considered, in various 
levels of approximation, the correlations introduced by 
the Coulomb interaction among charge carriers. It is well 
known that the physics of confined charge carriers is fun- 
damentally affected by the Coulomb interaction between 
charge carriers. As a consequence, a complete theoret- 
ical study of the physical properties of charge carriers 
in confining settings must take into account the full and 
non-trivial correlations that arise from the long-ranged 
Coulomb interaction. 

With this in sight, we present in this contribution a 
procedure to efficiently compute the matrix elements of 
the Coulomb interaction between charge carriers confined 
by a two dimensional parabolic potential and which are 
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FIG. 1: Bilayer of interplane distance 
at ri and r2, and relative angle 0. 



with charge carriers 
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spatially separated by a general interplane distance. This 
particular type of physical system is relevant in the study 
of graphene bilayers [5j, electron bilayers [9 indirect biex- 
citons [10] , bilayer quantum Hall systems pjj , self assem- 
bled quantum dots [l2[ , coupled quantum wells [H, Q , 
quantum dots formed by lateral fluctuations in the well 
[151 ], among many others. Having an efficient way to 
compute these elements opens the way for very relevant 
finite system calculations (i.e. 0(10) — O(10 2 ) particles), 
as opposed to the types of calculations that invoke the 
thermodynmic limit and that are usually encountered in 
some theoretical formalisms. Among the important finite 
system calculations, we may mention the Hartree-Fock, 
the Hartree-Fock-Bogoliubov (e.g. through the use of 
BCS-type wave functions) and the Random Phase Ap- 
proximation approximations [HI, [13] • 



In the case of the matrix elements of confined charge 
carriers moving in the same plane, the problem can be 
treated in an analytic form, yielding a g eneral expres- 
sion for the Coulomb matrix elements 18, 191. However, 
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when this formula is used for elements that involve states 
with high quantum numbers of either angular momenta 
or of radial excitations, convergence problems arise; it is 
then necessary to use a ppr opriate numerical methods to 
obtain reliable results [20, llj E2, [23[. We shall see in 
this paper that the computation of the matrix elements 
for the spatially separated case also involves the use of a 
special numerical procedure in order to accelerate the in- 
finite series that arises in the expansion of the expression 
for the Coulomb matrix elements. 

This paper is organized as follows: in Sec. H we review 
the Coulomb elements when the spatial separation d is 
equal to zero. In Sec. [Ill we give the expression for the 
Coulomb matrix elements calculated for the d ^ case, 
followed by Sec. IIIII where the method used to accelerate 
convergence is explained. Finally, in Sec. [iVlwe present 
some main results, and some conclusions are given in Sec. 
IVl Appendix I VII explains the notation and basis ordering 
used throughout the article. 



I. COULOMB MATRIX ELEMENT WITH d = 



Many authors have found various expressions for these 
elements [HI, [l9|, [26|. For example, using the two- 
dimensional Fourier transform for the Coulomb potential 
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it is possible to express the Coulomb matrix element in 
terms of finite sums in the form [l8[ 
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The eigenfunctions \i) = \ni,k) of a particle subject 
to either a two-dimensional harmonic potential [24] or a 
perpendicular magnetic field [25[ are given by 

W = = C nuU r^L^} (r 2 ) e^e" 1 *, (1) 

with m e {0}uN, k e Z, L^'(r) are the usual associated 
Laguerre polynomials and the normalization constant is 
given by C nu \ h \ = \Jni\j [tt (m + \k\)\]\ all lengths are 
scaled by the relevant length scale of the system i.e. the 
harmonic oscillator length or magnetic length, a choice 
which depends on the external potential being consid- 
ered. The corresponding energies are 



The Si ini i out expresses the conservation of the angular 
momentum l out = U + lj = h + h = hn- We have, 
for convenience, defined the quantities 1^ = 1^-/^1 = 
Ijh a ik = (m + m k + (\h\ + \l k \ - hk)/2) and A = 
((-irVm,!)(|Z,| + n^./mi + m)\(m - m t )\). 

Although expression (|6|) involves finite sums, it is inef- 
ficient for computing matrix elements that involve states 
of high angular momentum. The main difficulty arises 
from the computation of the large factorials. This ren- 
ders the calculation of the matrix elements computation- 
ally expensive. 

An alternative way of calculating the Coulomb matrix 
element is to write 
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for the Landau basis. 

Let (ij\V(d)\kl) = V(i,j,k,l,d) denote the Coulomb 
interaction matrix element for spatially separated charge 
carriers, with d the length of the separation (see Fig. 
[1]). Then, the matrix element between any two pairs of 
single-particle states for the d = case is written as 



(ij\V(0)\kl) = (ijh 
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It should be noted that is a 2D in-plane vector. 



By expanding the exponential inside the integral, replac- 
ing the result in (j4]) and integrating the angle and radial 
variables, we arrive at an explicit infinite series. The re- 
maining step of calculating the infinite series was tackled 
successfully in [27[ by making use of a series acceleration 
algorithm (different from the one used in this work). 



II. THE BILAYER CASE: COULOMB MATRIX 
ELEMENT WITH d ^ 

Let us now proceed to derive an expression for the 
general case d ^ 0. The Coulomb interaction in this case 
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reads 
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being the angle between the 2D in-plane vectors i*i and 
r 2 (Fig. [Tj), and the vector d being defined as d = dk. 
Note that ri, r 2 , and ri — r 2 are perpendicular to d. 
Manipulations analogous to those in the previous case, 
and use of Ll! 1 (r) = ^ =0 ((-l)-/m!)(|/| +n)!r-/((|/| + 
m)\(n — m)\) lead to the following expression for the 
Coulomb matrix elements 
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where 21 = \h\ + \lj\ + \lk\ + \h\, V = mi + rrij+mk + mi and 
iFi(a,6, z) is the confluent hypergeometric funtion (28[. 
This expression reduces in an equivalent form to (|6|) when 
d — » [23]. Although our expression involves an infinite 
sum, we can apply a numerical convergence procedure in 
order to make the computation of the matrix elements 
more tractable. To the best of our knowledge, there is 
no expression nor procedure for a practical calculation 
of the Coulomb matrix elements of spatially separated 
charge carriers in a two dimensional confining potential. 
As mentioned in the introduction, this result is of sig- 
nificant practical importance in finite system studies in 
various condensed matter systems. 

Incidentally, we would like to note first that it is possi- 
ble to reduce the matrix elements satisfying k — lj = h—h 
to only radial integrals of complete elliptic functions of 
the first kind K(k), as long as d > 0. Such a simplifi- 
cation is not possible when d = mainly because of the 
singular behaviour of the Coulomb interaction in such a 
case. The resulting integrals read 

(ij\V(d)\kl) = 

fiiinJout / dr i / dr2G(i,j,k,l,n,r2]d), (10) 
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Although the radial integrals have convergence problems 
for high values of the radial quantum numbers, due to 
the oscillatory nature of the Laguerre polynomials, the 
elements that can be computed agree with our results 
which make use of the series acceleration algorithm to be 
discussed presently. 



III. THE SERIES ACCELERATION 

We calculate the infinite series in the p index that ap- 
pears in (j9| using a u-type Levin transform to accelerate 
the convergence [29|. The basic idea is to construct an 
alternate series which converges much faster than the one 
we have at hand. Intuitively, this could be done if there 
was a way to simulate the asymptotic behaviour of the 
remainder R m of the series for large values of m. We can 
achieve this by constructing functions u m such that 
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c being of order unity. Having done this, we can then 
write 



(13) 



such that the coefficients fi m satisfy lim m ^ 00 /x m = c. By 
doing this, the asymptotic behaviour of the series will 
have been coded into the co m quantities. Now we must 
find a prescription for the coefficients fi m . This can be 
done by writing, for large m, the expression 



i=0 



Cii/Ji(m). 



(14) 



This expression is reasonable so long as the functions 
ipi(m) satisfy three conditions: first, we must have 
ipo(m) = 1 and Co — c Vm, so that ([T2]) holds; sec- 
ond, for i > we must have ^(m) — > when m — > 00; 
and third we must require that ^i+i(ra) = 0(^(m)) , so 
that, in taking the large m limit, all terms 2 / vanish 
at the same rate. The constants C{ are still unknown at 
this stage of the derivation. 

With these definitions, the m th partial sum s m of our 
series can be written as, for large m, 
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Here, s is the exact value of the series we want to com- 
pute. Finally, we truncate the sum, thus eliminating 
asymptotic terms: 

fc-i 

cr m = cr + uj m ^2 Cii/Ji(m). (16) 
i=0 

The corresponding change of notation from the exact s 
to the truncated a should be clear. Notice that the ex- 
pression ([T6|) tells us that the value a is such that all fc + 1 
equations, for m running from mtom + fc, hold simulta- 
neously. That is to say, the truncated value a depends on 
all terms which add up to a itself and, thus, it is as if it 
"feelt" the overall behaviour of the most important terms 
of the series. In this sense, the value calculated through 
the use of the u-type Levin transform can be thought of 
as an extrapolation of the final sum that uses informa- 
tion from only the first few terms of the complete series. 
The more terms we keep in the truncation, the better the 
extrapolation. 

In order to find the approximation a to the full sum, 
it is necessary to solve the fc + 1 equations for the fc + 1 
unknowns a and {Ci}. Hence, the resulting problem of 
calculating the series has been reduced to computing the 
solution of a fc + 1 dimensional matrix inversion problem. 
The larger the value of fc, the better the approximation 
to the exact value. The u-type Levin transform uses the 
choices ipi(m) = (m + /3)~i [29] and uo n = na n , which 
have proven to behave well for a large family of series [3(j • 
The parameter (3 is a real number that can be chosen to 
improve the rate of convergence. 

As an example, in our calculations, an element such as 

(15|VXd = 0.2)|23) « -0.0747592223 (17) 

required 1635 terms of the series to achieve convergence 
up to the tenth decimal place using the term-by-term 
sum, whereas it required less than 100 terms using the 
Levin transform, and was calculated in half the time us- 
ing a standard linear solver. It should be noted that the 
number of terms needed to achieve convergence will de- 
pend on the matrix element that is being calculated and 
that the calculation time will depend on the machine and 
the procedure used to solve the linear system that results 
from the Levin transform. 



IV. MAIN RESULTS 
A. Comparison between the d — and d ^ cases 

We first evaluate some particular matrix elements in 
order to get a feeling of how they vary as functions of the 
distance d. 

In Fig. ([2]), we show the behavior of the Coulomb ma- 
trix elements of the form (i, i|V(d)|i, i) as functions of 
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d 

FIG. 2: Elements of the form (i,i\V(d)\i,i) for: (Top) In 
the ordering of the harmonic oscillator basis. (Bottom) In 
the ordering of the Landau basis. The parentheses denote the 
quantum numbers (n, /). 



the separation between the planes, for both the oscilla- 
tor and Landau basis. For both basis sets, there is a 
general decaying behaviour. Clearly, when d ^> 1, all 
elements must approach zero. Refer to the appendix for 
the notation used for the basis sets. 

Although, there is a one-to-one correspondence be- 
tween both basis sets, there is a physical reason which 
justifies showing the calculations in both basis. The Lan- 
dau basis corresponds to the solution of a charged parti- 
cle moving in an external magnetic field perpendicular to 
the plane of motion. This system is infinitely degenerate 
with respect to the angular momentum quantum num- 
ber and, hence, its ordering is conceptually different to 
the finite number of degenerate energy levels of the two 
dimensional harmonic oscillator system. Thus, showing 
the calculation in each basis set is relevant for each of the 
two different types of in-plane confinement. 

Next, we construct a set of elements which involves 
pairs of states with fixed total angular momentum (i.e. 
using sets of two-particle states with some predetermined 
value of their total angular momentum). In order to do 
this, such states were randomly chosen for the cases I = 
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FIG. 3: Behaviour of the Coulomb average of randomly cho- 
sen pairs of states of total angular momentum / = (top) 
and I = 1 (bottom) as a function of d. This was done with 
the oscillator ordering. 



and 1 = 1. This is shown in Fig. ([3]). The general 
exponentially decaying behaviour observed in these plots 
seems to suggest that there is a scaling behaviour behind 
elements with predetermined total angular momentum. 

This last observation leads to an interesting behaviour 
which is worth emphasizing. Although all elements get 
suppressed in general as the distance is increased, they do 
not all approach zero at the same rate. Such a behaviour 
can be seen clearly in the top plot of Fig.©: one of 
the curves crosses two other curves, indicating that it is 
decaying faster than the other two. This tells us that, as 
d varies, the entries in the matrix representation of the 
Coulomb interaction vary in nontrivial ways with respect 
to each other. In the very large d limit i.e. when d ^> 1, 
we may obtain an expansion for the rate of decay 

d 1 

— V(iJ : k,l,d) « --pSikSji, (18) 

so that not all elements tend asymptotically to zero at 
the same rates, the rate depending on the indices of the 
states used in the matrix element. Hence it is possible, by 
changing the distance between the planes, for a system 
of charge carriers to exhibit different types of dynamical 
behaviour, as opposed to just interacting in a kind of 
screened Coulomb potential produced by the separation. 



In Fig. we show the most relevant elements for the dy- 
namics of Coulombian systems, namely (ij\V(d)\ij) and 
(ij\V(d)\ji), the so-called direct and exchange terms, us- 
ing the first energy shell of the Landau basis. Note how, 
as we increase the interplane distance, the direct terms 
get more suppressed, relative to the d = value, than 
the exchange ones. Also, the exchange elements decay 
faster than the direct ones as we change the index z, a be- 
haviour that is general for every d. This decay is expected 
because the exchange elements involve integrals of am- 
plitudes which can interfere destructively when summed 
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FIG. 4: Top: Direct Coulomb matrix element of the form 
V(i, 20, i, 20, d). Bottom: Exchange Coulomb matrix ele- 
ment of the form V(i, 20, 20, i, d). Both plots were done with 
the index i running through the first 60 states of the first 
Landau level 



over, as opposed to the direct elements which involve 
probability densities which are always positive. Finally, 
the maximum value of the direct term is actually shifted 
from its diagonal value V(i,i,i,i,d), and this maximum 
gets further shifted as d increases. However, the maxi- 
mum for the exchange plot remains at the diagonal point, 
although this is actually again the same diagonal (and, 
thus, direct) matrix element V(i,i,i,i,d). 



B. Amplitude transitions and diagonaliation of a 
single indirect exciton 

Since the Coulomb matrix elements connect two pairs 
of two-particle states through the Coulomb interaction, 
we can use them to weigh how probable it is for a partic- 
ular transition to occur as a function of the distance d. In 
order to study a simple transition, we computed the ma- 
trix elements for a given pair of initial states and several 
pairs of final states, all pairs having zero total angular 
momentum. This is shown in Fig. ([5]). Each continuous 
line is the absolute value of the transition amplitude for 
a given distance d and each point represents a final state, 



6 




Transition 




Transition 

FIG. 5: Absolute value of the transition amplitudes between 
harmonic oscillator states of zero total angular momentum. 
Each continuous line has a fixed d value. As d increases, the 
lines get suppressed. The initial state corresponds to two 
particles in |1)|1) (top) and |4)|6) (bottom). The final two- 
particle states are ordered horizontally according to the list: 
|1)|1>, |1}|5), |2)|3>, |3)|2), |5)|1), |2}|9>, |3}|8>, |4)|6>, |5)|5), 
|6)|4>, |8)|3>, |9)|2>, |7)|10>, |8)|9>, |9)|8>. 



as explained in the caption of the figure. The succession 
of lines go from d = through d = 1. As d increases, 
all points (i.e. all considered transitions) generally get 
suppressed. 



the unequal rates of change as d is varied only this time 
it is evidenced with the eigenenergies of the exciton. As 
d increases, the spectrum tends to the harmonic oscil- 
lator energies for two particles with total zero angular 
momentum, as expected. 
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FIG. 6: First few exciton energies Ei for zero total angular 
momentum. Each continuous line represents a fixed value of d. 
The inset shows the binding energy, defined as Eh — E\—E osc 
(i.e. the portion of the energy of the bound electron-hole pair 
which is due solely to the Coulomb interaction) , as a function 
of the interplane distance d 



Once again, we note that the rate at which each coeffi- 
cient gets suppressed is not the same for all, the diag- 
onal transition being typically the most significant one, 
as further calculations of other similar cases have shown. 
In fact, when d is of order 1 (i.e. when the distance is 
of the order of the oscillator length), almost all of the 
coefficients are negligible, except for the diagonal ones. 
This reasonably shows that, when the planes are sepa- 
rated beyond the characteristic distance of the system, 
all transitions get suppressed, albeit at unequal rates. 

Finally, as a more physical application, we have consid- 
ered the simple case of an electron and a hole spatially 
separated in an effective mass hamiltonian with a har- 
monic potential. The eigenstates of this system are those 
of a single spatially indirect exciton. The hamiltonian, 
in dimensionless units, reads 

H = Y,i4 e) 4ei + 4 h) h\hi} -I3j2(ij\\kl)elhihje k . 

1 ij kl 

(19) 

In this expression, we have made the convenient defini- 
tion (ij\\kl) = (ij\V(d)\kl). The (3 parameter denotes 
the ratio between the Coulomb and harmonic oscillator 
energies, and the ef* and e-^ are the dimensionless har- 
monic oscillator energies of the electron and the hole, 
respectively. The diagonalization for the I = block of 
the hamiltonian is shown in Fig.©. As before we notice 



C. Convergence sum rule 

As a consistency check of the above results, we com- 
puted an identity involving the Coulomb matrix el- 
ements, reminiscent of the sum rules encountered in 
atomic physics. Namely, it can be checked that, pro- 
vided li — lj = Ik — k, we can integrate out the angular 
integrals exactly to obtain a relation of the form 

{ij\V\d)\ij)=J2\(ij\V(d)\kl)\ 2 = (20) 

kl 

noo roo 

$iin,iou t / dr i / dr 2 F(i,j,k,l,n,r2]d), 
Jo Jo 

with 

F (i,j,k,l,r 1 ,r 2 ]d) = Cn^Cn^C^C^ 
Xr \li\Ml k \+^ r \h\Mh\ + l e -{rl+rl) 

< 2 " 2 (21, 
^f{<P + (n - r 2 ) 2 )(d 2 + (n + r 2 ) 2 ) 
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FIG. 7: Convergence of the left hand side of Eq. (|20|) as a 
function of the number of oscillator states \kl) used in the 
summation, starting from the two particle states which are 
closer to the state \ij) (the i in horizontal axis is an index that 
denotes the i th element of the two-particle oscillator basis of 
the form \kl)). In this particular example, the state \ij) is 
the ground state. The horizontal dashed lines are the exact 
values calculated with the right hand side of Eq. (|20[) . 

Note that the convergence of the curve corresponding 
to d = 0.2 is slower than that of the other curves. This 
has a simple and clear physical meaning. As the spatial 
separation between the planes becomes smaller, it will be 
easier for the Coulomb interaction to induce transitions 
starting from two particles in the oscillator ground state 
to two particles in two other higher oscillator states. This 
means that all amplitudes that we add in the sum rule 
become more significant as d becomes smaller and, thus, 
the sum will require more and more terms which involve 
higher oscillator states. Notwithstanding this, there is 
clear convergence. 

The radial integrals can be evaluated numerically by 
various methods, which makes the computation of the 
right hand side of ([20]) completely independent from our 
main numerical approach, which is in turn used to com- 
pute the left hand side. We found that the relation ([20]) 
was very well satisfied, as can be seen in Fig. ([?)). 



V. CONCLUSIONS 

In the work presented here, we have shown the pro- 
cedure and results of a numerical calculation of the 
Coulomb matrix elements of spatially separated charge 
carriers confined either by a two dimensional parabolic 
quantum dot or through the use of a magnetic field. We 
have found that, in a very important way, each matrix 
element is affected differently as d is varied. This could 
lead to different and interesting physical regimes of a sys- 
tem subject to interparticle Coulomb interaction. 

Although the procedure we have implemented above 



. . . (2,-3) (2,-2) (2,-1) (2,0) (1,1) (0,2) 

|6 3 ) |5 3 ) |4 3 ) |3 3 ) |2 3 ) |1 3 ) 

. . . (1,-3) (1,-2) (1,-1) (1,0) (0,1) 

|5 2 ) |4 2 ) |3 2 ) |2 2 ) |1 2 ) 

. . . (0,-3) (0,-2) (0,-1) (0,0) 

|4i) |3i> |2 X ) 

FIG. 8: Ordering of the Landau basis. Each line is an energy 
shell which has infinite many elements. The notation \Nu) 
refers to the N th state of the // Landau level. 




FIG. 9: Ordering of the harmonic oscillator basis. Each line 
is an energy shell which has finite many elements. 



has worked quite well, it should be noted that, as with 
any numerical algorithm, care must be taken in evalu- 
ating elements with too high a difference in quantum 
numbers of the involved states e.g. high differences be- 
tween the U and Ik or kj and k\. These cases can become 
pathological because of the highly oscillatory nature of 
the associated wavefunctions. 

We would like to emphasize the relevance of the pro- 
cedure presented here. With such elements, theoretical 
studies can be carried out that are comparable with ex- 
perimental systems. Furthermore, this study is relevant 
for describing and understanding many-body effects such 
as collective behaviour and quantum correlations. 

Work that involves the in-depth study of spatially 
indirect excitons within a finite system framework using 
Hartre-Fock and BCS type approximations is already 
on the way. For such a study, the matrix elements 
computed in the present paper will be critically needed. 
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VI. APPENDIX: BASIS ORDERING 

Since we will be refering to two different basis sets, 
we will specify the ordering used throughout the pa- 
per. Both basis sets are organized in shells according 
to their respective energies given in Eq.([2]) and Eq.(j3]). 



For the Landau basis, the logic of the numbering can be 
extracted from Fig.®; and that of the oscillator basis, 
from Fig.©. The parentheses in these figures follow the 
notation (n, I) i.e. the radial and the angular momen- 
tum quantum numbers associated to the corresponding 
wavefunctions shown in Eq.([T]) of the charge carriers. 
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